subroutine triangulation (PointM, Nonvalue,NODE_CHANGE, n_regrid, m_regrid, resolx, resoly, xmin, xmax, ymin, ymax)

implicit none

integer  n_regrid,m_regrid, cv, tv, NODE, Coordiinf, Coordjinf, rowmin, i, k, j, bestconfig
parameter (NODE=6100)
parameter (cv = 2)
parameter (tv = 2)
real  PointM(3), coordi, coordj, Nonvalue, resolx, resoly, C, B, A, DELTA
real  NODE_CHANGE (,)

real  coordireal, coordjreal, xmin, xmax, ymin, ymax, dist
real, dimension (4,4)  distance
real  D(6), E(2), Projection(3), Triangle(2), Sumtriangle(2)
real  PointA(3), PointB(3), PointC(3), vectorAB(3), vectorAC(3), vectorN(3)



coordi = PointM(1)
coordj = PointM(2)

Coordiinf = coordi
Coordjinf = coordj

distance(1,1)=xmin+(Coordiinf-1)Resolx
distance(1,2)=ymin+(Coordjinf-1)Resoly
distance(1,3)=NODE_CHANGE(Coordiinf,Coordjinf)
    
distance(2,1)=xmin+(Coordiinf)Resolx
distance(2,2)=ymin+(Coordjinf-1)Resoly
distance(2,3)=NODE_CHANGE(Coordiinf,Coordjinf)
    
distance(3,1)=xmin+(Coordiinf-1)Resolx
distance(3,2)=ymin+(Coordjinf)Resoly
distance(3,3)=NODE_CHANGE(Coordiinf,Coordjinf)
    
distance(4,1)=xmin+(Coordiinf)Resolx
distance(4,2)=ymin+(Coordjinf)Resoly
distance(4,3)=NODE_CHANGE(Coordiinf,Coordjinf)

if(distance(1,3)-1Nonvalue.or.distance(2,3)-1Nonvalue.or.distance(3,3)-1Nonvalue.or.distance(4,3)-1Nonvalue) then;
PointM(3)=Nonvalue
else;
do i=1,4
distance(i,4)=SQRT((Coordi-distance(i,1))2+(Coordj-distance(i,2))2)
end do


    dist=distance(1,4)
    rowmin=1
	do i=2,4
    	if (distance(i,4)dist) then;
        	dist=distance(i,4)
            rowmin=i
        end if
    end do
    
D(1)=SQRT((distance(3,1)-distance(2,1))2+(distance(3,2)-distance(2,2))2+(distance(3,3)-distance(2,3))2)
D(2)=SQRT((distance(2,1)-distance(1,1))2+(distance(2,2)-distance(1,2))2+(distance(2,3)-distance(1,3))2)
D(3)=SQRT((distance(1,1)-distance(3,1))2+(distance(1,2)-distance(3,2))2+(distance(1,3)-distance(3,3))2)
D(4)=SQRT((distance(4,1)-distance(3,1))2+(distance(4,2)-distance(3,2))2+(distance(4,3)-distance(3,3))2)
D(5)=SQRT((distance(4,1)-distance(2,1))2+(distance(4,2)-distance(2,2))2+(distance(4,3)-distance(2,3))2)
D(6)=SQRT((distance(4,1)-distance(1,1))2+(distance(4,2)-distance(1,2))2+(distance(4,3)-distance(1,3))2)
    
    
do k=1,2        
		do j=1,2
    	E(j)=((D(5k-4))2+(D(2j-k+2))2-(D(6j+7k-4jk-7))2)(2D(5k-4))
		do i=1,3
        Projection(i)=distance(10-5j-6k+4kj,i)
        Projection(i)=Projection(i)+((distance(5j+6k-4kj-5,i)-distance(10-5j-6k+4kj,i))E(j))D(5k-4)	
		end do
		
		E(j)=(Projection(1)-distance(2k+j-2,1))2+(Projection(2)-distance(2k+j-2,1))2
        E(j)=E(j)+(Projection(3)-distance(2k+j-2,1))2
        E(j)=SQRT(E(j))
        
		Triangle(j)=E(j)D(5k-4)2
		end do
		Sumtriangle(k)=Triangle(1)+Triangle(2)
	end do

	bestconfig=1
	if (Sumtriangle(2)Sumtriangle(1)) then;
  		bestconfig=2
	end if

	if (bestconfig==1) then;
    	    do i=1,3
			PointB(i)= distance(3,i)
            PointC(i)= distance(2,i)
            end do
            rowmin=1
            if (distance(1,4)distance(4,4)) then;
            rowmin=4
            end if
    else;
   			do i=1,3
            PointB(i)= distance(1,i)
            PointC(i)= distance(4,i)
            end do
            rowmin=2
            if (distance(2,4)distance(3,4)) then;
            rowmin=3
            end if
    end if
    do i=1,3
    PointA(i)= distance(rowmin,i)
    end do
    
	do i=1,3
		vectorAB(i)=PointB(i)-PointA(i)
	end do

    do i=1,3
      	vectorAC(i)=PointC(i)-PointA(i)
    end do

	vectorN(1)=vectorAB(2)vectorAC(3)-vectorAB(3)vectorAC(2)
    vectorN(2)=vectorAB(3)vectorAC(1)-vectorAB(1)vectorAC(3)
    vectorN(3)=vectorAB(1)vectorAC(2)-vectorAB(2)vectorAC(1)
    
	A=vectorN(1)
    B=vectorN(2)
    C=vectorN(3)
    Delta=-PointA(1)vectorN(1)-PointA(2)vectorN(2)-PointA(3)vectorN(3)

    if (C==0) then;
    PointM(3)=Nonvalue
	else;
    PointM(3)=(APointM(1)+BPointM(2)+Delta)(-C)
	end if

   
end if
    
return
end subroutine triangulation